clear
clc

E = 1; % энергия сигналов

numFig = 0; % счетчик графиков

trel = struct('numInputSymbols', 2, 'numOutputSymbols', 4, 'numStates', 2,...
              'nextStates', [0 1; 0 1], 'outputs', [-1 1; -1 1]);

m = log2(trel.numInputSymbols);

k = 10; % число бит информации

% Parameters of the channel generator
F        = [0.1,0.5];  %[f0,f1,...]
mCount   = 2;         % Mirrors amount
timeStep = 0.01;      % Time to send one bit
testMsg  = [1,1,-1];  % Test signal

%% Моделирование

nErrMax = 990; % число ошибок, при котором завершается моделирование при
               % данном значении отношения сигнал/шум. Происходит переход
               % к следующему значению отношения сигнал/шум

listSNRdB     = [18 20 22 24 25 26 27 28];
listSNRdB_bit = listSNRdB - 10*log10(k);
Pe_bit        = zeros( size(listSNRdB) );

%% Моделирование
j = 1; % счетчик отношений сигнал к шуму
for SNRdB = listSNRdB
    
    
    % SNR = E/No
    % SNR - signal to noise ratio (отношение сигнал/шум)
    SNR = 10^(SNRdB/10); % перевод из децибел
    
    noise  = E/SNR;        % шум
    DNoise = noise/2;      % дисперсия шума, для БГШ D = NO/2
    sigma  = sqrt(DNoise); % среднее квадратичное отклонение
    
    
    tStart = 0;  % Start time
    
    l = 0; % счетчик испытаний
    
    nBitErr = 0; % счетчик битовых ошибок
    
    fprintf('SNRdB = %f\n', SNRdB);
    
    while nBitErr < nErrMax
        % псевдослучайные данные
        msg = [1; randi([0 1], [k-1, 1])]; % обеспечить начальное состояние решетки детектора
        
        % двоичная модуляция
        points = msg.*2 - 1;
        
        % Time interval definition
        tStop = tStart + (k-1 + length(testMsg))*timeStep;
        t = tStart : timeStep : tStop;
        
        
        % Mu coeffs generetion
        MUcoef = raylrnd(ones(mCount,size(t,2)));
        %MUcoef = ones(k,size(T,2))
        MU = (MUcoef.*0.0001 + 1) .* (cos(2*pi*F'*t).*0.8 + 1);

        
        % добавление тестовых битов
        points = cat(2,points',testMsg);
        
        Channel = zeros(1,length(points));
        
        % Channel modeling
        Channel(1,1)=points(1,1)*MU(1,1);
        for i = 2:length(points)
            Channel(1,i) = points(1,i)   * MU(1,i) + ...
                           points(1,i-1) * MU(2,i);
        end
        
        % оценка параметров канала
        A = [testMsg(1,2) testMsg(1,1)
             testMsg(1,3) testMsg(1,2)];
        f = [Channel(1,k+2); Channel(1,k+3)];
        MuEst = A\f; % Estimated channel coeffs mu0 mu1
        
        % Time shift for the next iteration
        tStart = tStop + timeStep;
        
        % удаление тестовых битов
        Channel = Channel(1,1:k);
        
        % добавить шума
        for i = 1 : k
            Channel(i) = Channel(i) + sigma*randn;
        end
        
        % обновить решетку кода
        trel.outputs = [ -MuEst(1)-MuEst(2) MuEst(1)-MuEst(2)
                         -MuEst(1)+MuEst(2) MuEst(1)+MuEst(2) ];
        
        decoded = myvitdec(Channel, trel);
        
        if ~isequal(msg, decoded)
            nBitErr   = nBitErr + biterr(msg, decoded);
            Pe_bit(j) = nBitErr / (l * k); % уточнить вероятность ошибки
            
            % уточнить график
            semilogy(listSNRdB_bit, Pe_bit, 'k-o');
            title('P_{e}^{(bit)}(E/N_{o})_{(bit)}')
            xlabel('(E/N_{o})_{(bit)}, dB')
            ylabel('P_{e}^{(bit)}')
            grid on
            drawnow
        end
        
        l = l+1;
    end
    
    j = j+1;
end

semilogy(listSNRdB_bit, Pe_bit, 'k-o');
title('P_{e}^{(bit)}(E/N_{o})_{(bit)}')
xlabel('(E/N_{o})_{(bit)}, dB')
ylabel('P_{e}^{(bit)}')
grid on
